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Abstract 

In this paper we present a network model to study the impact of spatial distribu- 
tion of constituents, coupling between them and diffusive processes in the context 
of biological situations. The model is in terms of network of mobile elements whose 
internal dynamics is given by differential equations exhibiting switching and/or os- 
cillatory behaviour. To make the model more consistent with the underlying biolog- 
ical phenomena we incorporate properties like growth and decay into the network. 
Such a model exhibits a plethora of attributes which are interesting from both the 
network theory perspective as well as from the point of view of bio-chemistry and 
biology. 

We characterise this network by calculating the usual network measures like net- 
work efficiency, entropy growth, vertex degree distribution, geodesic lengths, cen- 
trality as well as fractal dimensions and generalised entropy. It is seen that the 
model can demonstrate the features of both scale free networks as well as small 
worlds network in different parameter domains. The formation of coherent spatio- 
temporal patterns is another feature of such networks which makes them appealing 
for understanding broad qualitative aspects of problems like cell differentiation (e.g. 
morphogenesis) and synchronization (e.g. quorum sensing mechanisms). 

One of the key features of any biological system is its response to external attacks. 
The response of the network to various attack strategies (isolated and multiple) is 
also studied. 
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1 Introduction 



Biological systems exhibit tremendous complexity, versatility and robustness 
in their responses to their environment. Though molecular biology has had 
spectacular success it is becoming clearer that simple enumeration of genes, 
proteins and metabolites is not sufficient to understand this complex be- 
haviour. In recent time it has been shown that [Lee et.al,Milo et.al] to cap- 
ture system level behaviour, it is instructive to search for patterns common 
to complex systems and networks in general. For example, even if we could 
calculate the physio-chemical properties of a protein from its structure, to ob- 
tain the physiological behaviour from these would be a daunting task without 
the knowledge about the organisation of the network in which the proteins 
operate. Networks have been discovered to be important at various levels in 
biological systems [Hartwell et.al.]. These networks of interacting elements are 
intrinsically dynamic and describe how the system changes in space and time 
in response to the external stimuli, grows and reproduces, differentiates and 
dies. Therefore there is a pressing need in biology to develop an alternative 
language and methodology to deal with such situations. 

"Perhaps a proper understanding of the complex regulatory networks making 
up cellular systems like the cell cycle will require a shift from common sense 
thinking. We might need to move into a strange more abstract world, more 
readily analyzable in terms of mathematics than our present imaginings of 
cells operating as a microcosm of our everyday world" [Nurse P.] 

"The best test of our understanding of cells will be to make quantitative pre- 
dictions about their behaviour and test them. This will require detailed sim- 
ulations of the biochemical processes taking place within [cells]. We need to 
develop simplifying, higher-level models and find general principles that will 
allow us to grasp and manipulate the functions of [biochemical networks]. " 
[Hartwell et.al.] 

This modelling of processes in cellular and molecular biology is based upon the 
assumption that interactions between molecular components can be approxi- 
mated by a network of biochemical reactions in an ideal macroscopic reactor. 
Then following the standard methods of chemical reaction kinetics, ordinary 
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differential equations are obtained. This modeling approach has been applied 
to many systems. 

However, there are a few considerations which are omitted in above approxima- 
tion. Firstly, though a few studies consider the spatial aspects of the problem 
(e.g. [Eldar et. al (2002)] on Drosophila melanogaster), most of the time they 
are simply neglected. However, it is known that the cell is not a homogenous 
well stirred reactor, but instead is a highly compartmentalized and heteroge- 
nous environment where phenomena like crowding or channeling are present 
[Ellis E.(2001)]. In such a situation, the spatial aspects may well play a key 
role and may need to be incorporated. Such spatial models are particularly 
relevant for biological networks where the concentration and growth of sig- 
nalling factors decays with the distance. Thus the probability of establishing 
a connection with far off nodes is much smaller. Thus the nodes in such net- 
works are more likely to be directly connected to the "nearby nodes" than in 
the far off nodes. The spatial inhomogeneity observed in real life networks is 
also indicative of the fact that the values in the network are not uniformly dis- 
tributed but instead clustered together at various points in space. There may 
also be bounds and limits in spatial model which forbid or allow connections 
depending on a cut-off distance. 

Another important consideration is that of mobility. In most of the studies 
the elements are considered to be fixed in space and the structure of coupling 
is fixed in time i.e. two elements that are initially coupled are coupled for- 
ever and they remain at the same point in space. Some models do consider 
time dependent coupling coefficients [Shibata T. and Kaneko K.] though spa- 
tial positions are fixed. However, in real situations the elements move due to 
diffusion or chemical gradients [Painter K.J and Hillen T.(2002)] and the local 
interactions determine the connectivity of the element. 

There is another way in which such models make restrictive assumptions in 
context of biological systems. The number of nodes and degrees of freedom are 
fixed initially and remain fixed. However most of the biological systems show 
growth, reproduction and death. For example a cell multiplies and transfers 
its attributes to its offspring and dies. There may be more than one pathway 
to produce a protein and also more than one way it is degraded. The biological 
elements also differentiate themselves to allow flexibility with respect to func- 
tion [Huang S. and Ingber D.E. (2000)]. This differentiation and growth leads 
to degrees of freedom no longer remaining fixed for a node and diversification 
of the nature of the elements themselves. 

Lastly, most of the models assume either a boolean function for evolution 
[Shmulevich I. et.al. (2002)] or an arbitrary dynamical rule like dynamical 
maps(e.g. logistic, circle or Cellular Automata update rules) or some conve- 
nient rule like piecewise linear differential equations [de Jong H. et. al.(2003)]. 
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However, such a dynamical rule should have some correspondence(e.g. same 
kind of non-linearity, continuity properties) to the underlying biological con- 
sideration so that it can reproduce behaviour more realistically. Ideally this 
function should be inferred from experimental data or the detailed knowledge 
of underlying kinetics. However, as an approximation, at least some motif bi- 
ological evolution rules which are found in a number of contexts should be 
used to model the broad qualitative features with a hope that they fall in the 
same class of systems as the one being modelled and reproduce results with a 
greater degree of faithfulness to reality. 

Thus in light of the above discussion, there is a need to generalise the existing 
network models of the biological systems. Here we attempt to relax all of the 
assumptions made above and study the behaviour of such models. 

The paper is organised as follows: In the next section we set up the network 
and elaborate on its features. We also discuss the internal dynamics of the 
nodes and various choices one can make. The section on analysis studies the 
mathematical properties of the network and compares them with numerical 
results. After that we consider various strategies of attack on the network to 
determine its robustness. Finally we conclude with discussion and interpreta- 
tion of results. 



2 Materials and methods 

In this section we set up the network as follows: 
2. 1 Defining Network Behaviour 

A network consists of a set of connected nodes. Each node has an internal 
state and an intra-dynamics which evolves the state according to a dynamical 
system. The current state of the node is also dependent on the state of all 
its neighbours by coupling. In the network that we set up, the number of 
neighbours a node is coupled to is dependent on the dynamics of the network 
(inter-dynamics). This is in contrast to the lattice models where the number of 
neighbours is fixed a-priori. In such lattice models, one either couples to a fixed 
number of nearest neighbours or has a global coupling to all the other nodes 
in the network (mean- field models). We feel that in a number of biological 
contexts all the above assumptions are a little too restricted. During the course 
of time evolution, a given sub part of the system usually interacts with a 
variable number of other elements [Shibata T. and Kaneko K.]. Furthermore 
the network that we set up is an evolving one, where new nodes get added and 
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redundant nodes get removed. Therefore, the number of neighbours should be 
dynamically changing. 

Evolving networks have been studied in recent times in great detail 
[Erdos et al(1961); Barabasi(2001); Watts(1998); Newman]. Typically the net- 
works evolve by adding nodes to a random initial configuration (e.g. Erdos- 
Renyi network model). The nodes are usally added by using different heuris- 
tic prescriptions. In particular, scale free networks are obtained by adding 
new nodes based on the vertex degree distribution. These networks have been 
subject of much study and debate [Albert, R. and Barabasi, A.-L. (2002). ] 
and model many biological contexts, for example phenomena involving food 
webs, biochemical interactions, protein-protein interaction maps for viruses, 
prokaryotes and eukaryotes C.elegans etc are described by scale free networks. 
A different prescription (e.g Watt's Beta model) results in a small worlds net- 
work which is found to be useful in modeling co-operative behavior in bi- 
ological systems [Erdos et al(1961); Barabasi(2001); Watts(1998); Newman]. 
However, due to finite size and limited data available in most realistic sit- 
uations, it is difficult to establish scale free behavior of a network unambigu- 
ously [Albert, R. and Barabasi, A.-L. (2002). ]. Secondly, the scale free nature 
is often inherent in the prescription to grow a network. 

We take a different approach in the present work. Rather than starting from a 
random network, we start with a network of few elements and then grow the 
network by state-splitting [Lind D. and Marcus B. (1995)]. Such a growth law 
does not have a built in prescription for power law behaviour. Thus any scale 
free nature found is truly an emergent property rather than a consequence of 
a heuristic for network growth. 

Diffusion plays a major role in most of the biological processes. Thus any model 
which describes a process involving spatial motion or distribution of biological 
elements should take into account this important feature. The existing network 
models try to incorporate this by varying the connections between the elements 
dynamically [Shibata T. and Kaneko K.]. However the assumptions that are 
made in these models are too restrictive in many situations. We describe these 
assumptions and their possible relaxations below. 

(1) State is mostly considered to be discrete. For example most com- 
monly used network for modelling is a boolean network in which internal 
state of the nodes is a discrete (0 or 1) state. This kind of model is, 
at best, an approximation for modelling qualitative coarse grained be- 
haviour of the system. Some of the limitations of a boolean approach are 
the following : firstly since the internal state of an elemnent is evolving 
in a discrete manner, it becomes difificult to model situations which in- 
volve elements with different times scale being updated asynchronously. 
Secondly, we can not have a coupling between elements which is weighted 
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by the continous value of the internal state. Such a coupling is desirable 
particularly in bio-chemical reaction networks Thirdly, as a dynamical 
system, the nodes can have more than one attractor i.e. more than one 
possible final state depending on initial values. To model such situations 
with boolean networks is cumbersome. On the other hands, using differ- 
ential equation has advantages of preserving " quantitativity and causal- 
ity" [De Hoon M. et. al.(2003)] In more realistic models, the nodes have 
a number of states possible. We assume the state variable describing the 
internal state to be a continuous variable thus allowing more quantitative 
approach. 

(2) The internal state of the node does not affect the interactions in 
the networks. In realistic situations this assumption is restrictive. This 
assumption is actually inherent in the nature of the dynamical systems 
models where couplings are usually considered as fixed. There are models 
which consider variable topologies 

[Erdos et al(1961); Barabasi(2001); Watts(1998); Newman]. Even in such 
models the nodes are usually kept fixed and the connections that they 
make are dynamically changing. However, in most of the situations the 
elements move (either along a signal gradient or in a diffusive manner) in 
a very heterogenous medium. Most of the interactions are also dependent 
on the actual concentrations of active features of these elements. These 
active features depend on the internal state of the node and change with 
time. Thus the number of other nodes a given element would interact with 
depends on the current values of its internal state as well as the current 
values of internal states of the other elements. Therefore there is a need 
for introducing an explicit coupling between the internal states and the 
strength of interaction. We address this issue as follows: The nodes are 
allowed to move according to either a diffusive motion or a gradient or 
both. The interaction that takes place depends upon the number of other 
nodes in its "sphere of influence" and the internal states of the interact- 
ing elements have an effect on each other. The mobility of the element is 
dependent on its internal state and the region it finds itself in. Thus we 
have an explicit coupling between the number of neighbours that a node 
has and its internal state. 

(3) The interaction range is independent of time. However, for bio- 
logical systems one can not categorise them as having a fixed interaction 
range but should have coupling terms which are changing with time. This 
is because such systems usually have a number of spatio-temporal scales 
with different elements interacting at different scales. This is modelled 
in the current system by using a coupling strength that is dependent on 
the distance between the interacting elements. This dependence can in- 
corporated either by using a threshold value or by using an interaction 
strength which decays with distance. Thus as an element moves around, 
its interaction range with other elements is dynamically changing. 

(4) The structure of interactions does not change with time i.e. if two 



6 



elements are coupled then they have no mechanism to decouple and two 
decoupled elements do not couple. We allow for such a mechanism. If an 
element is outside the sphere of influence its interaction diminishes and if 
it moves within sphere of influence of other element it gets coupled. Thus 
the overall network dynamics encompasses many coupling and decoupling 
events and overall evolution of network travels through many structures 
of interactions. 

(5) Though there is a growth law (fixed by network evolution heuris- 
tic) there is usually no decay law. As mentioned earlier elements 
can grow, reproduce and die in many situations of interest. The de- 
cay law, modelling the removal or death of elements in biology plays 
an important role in determining behaviour. We incorporate these ef- 
fects as follows: the growth is modelled by state splitting of the network 
[Lind D. and Marcus B. (1995)]. Such a mechanism is natural for these 
contexts. There are two kinds of state splitting that are considered. In- 
splitting in which the incoming degree of a node is split between the 
node and its offspring and the outgoing degree is shared by both and 
the out-splitting in which outgoing degree is split and incoming connec- 
tions are shared. This is different from self-organised scale free network 
in which nodes get all the degrees inherited. The decay is modelled as 
following: once a node is out of the interaction range of all other nodes, 
it decouples and is marked as dying. In such a state it would evolve as 
an independent dynamical system and reach a attractor of its evolution 
state and stay at that state without contributiong to the dynamics on the 
network. However, since other nodes are also mobile it can happen that 
another node with connections to the network comes within the sphere 
of influence of this node and get coupled to it. Then the node is "healed" 
and starts contributing to the network dynamics again. The nodes which 
do not heal within a specified time are considered dead and are removed 
from the system. Such a removal does not affect the dynamics of the rest 
of the system and is analogous to how biological systems clear dead and 
redundant parts. 

(6) The elements are considered identical. We incorporate diversity and 
differentiation of the elements in our system. Firstly, the network is made 
up of a number of different types of elements which are representative of 
various kinds of dynamics actually observed in biological networks (see 
next section). 

Thus our network is described by the following: We have a set of coupled nodes. 
The number of nodes each element couples to is determined by its state and 
its sphere of influence. Every node also has a variable defining its position 
in space. This variable depends on the current state of the element, current 
position of the element, current position of its neighbours and current state of 
the neighbours. Thus the interaction structure of the network is dynamically 
dependent on the states, positions and spheres of influence of the nodes. 
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The network grows by reproduction of its nodes. To model this growth we use 
state splitting and to model decay we use the following rule: if two nodes are 
greater than a certain distance apart then they are set to be decoupled. A 
node which has been decoupled from every node is a "dying" node. However, 
due to dynamic interactions, such a node may still be revived ("healed"). Any 
dying node which does not heal for a specified period is considered "dead" 
and removed from the network. 

Thus the network is described by the following set of equations 
Calculation of the state of a node. 



Xf 1 = X$ + (1/Nl 3 ) £ (e m (i)f [ Hx\)) ~ {l/Nl 3 ) £ (e out (i)fM(X\)) (1) 

i=i i=i 



Calculation of the current position of a node. 
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Calculation of the number of in-neighbours and out-neighbours of a node. 
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where 



X* — > The state vector for a node j at time t 

r* — > The position of a node j at time t 

N[ — > The number of in neighbours of a node j at time t 

jV| — > The number of out neighbours of a node j at time t 

e m (i) — >■ The in coupling coefficient for node i 

Cout(i) — ¥ The out coupling coefficient for node i 
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c(R) — > is a coefficient depending on radius of influence R with value between 
and 1. 

e(\ri- rj \)=0 if\n-r 3 \ > R 

= 1 otherwise (5) 

2.2 Defining Node Behaviour 

Every node has an internal state and an internal dynamics. For the internal 
dynamics we choose some well known motifs which are frequently seen in dy- 
namics of signalling and regulatory pathways in the cells [Tyson (2003)]. These 
motifs, which behave as switches and oscillators, serve as building blocks for 
the bigger, more complex networks. For example, consider a typical situation 

The proteins that modulate Cdc2 activity are themselves modulated by 
Cdc2:Cdcl3, through a set of feedback loops 

(1) Ruml inhibits Cdc2:Cdcl3, but Cdc2:Cdcl3 phosphorylates Ruml, thereby 
targeting Ruml for degradation. 

(2) Ste9:APC labels Cdcl3 for degradation, but Cdc2:Cdcl3 can phosphory- 
late Ste9, thereby downregulating its activity and targeting it for degra- 
dation. 

(3) Weel phosphorylates and inactivates Cdc2:Cdcl3, but, at the same time, 
Cdc2:Cdcl3 is trying to phosphorylate and inactivate Weel. 

(4) Cdc25 takes the inactivating phosphate group off PCdc2: Cdcl3, and 
Cdc2:Cdcl3 returns the favor by phosphorylating and thereby activating 
Cdc25. 

(5) Slpl:APC, which also labels Cdcl3 for degradation, is itself activated by 
Cdc2:Cdcl3 by an indirect pathway. 

The first three feedback loops are examples of mutual antagonism. Under 
appropriate conditions, the antagonists cannot coexist, i.e. the feedback loop 
works like a toggle switch. Either Cdc2:Cdcl3 has the upper hand and its 
antagonist (Ruml or Ste9 or Weel) is suppressed, or vice versa. The fourth 
interaction is a positive feedback loop: Cdc2 and Cdc25 activate each other in 
a mutually amplifying fashion. The last interaction is a time-delayed negative 
feedback loop, which, under appropriate conditions, can generate oscillations 
(as Cdc2:Cdcl3 concentration rises, it turns on Slpl, which targets Cdcl3 for 
degradation, causing Cdc2:Cdcl3 concentration to fall, and Slpl to turn off) 

To study the dynamical consequences of these feedback loops, one must formu- 
late these interactions as a precise molecular mechanism, convert the mech- 
anism into a set of nonlinear ordinary differential equations, and study the 
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solutions of the differential equations by numerical simulation. 



We can set up a network by randomly selecting these elements and connecting 
them. However, for the present study we choose four elements behaving as 
switches and two behaving as oscillators. 

The equations for the switches and oscillators that we use are as follows 
Perfectly Adapted Switch 

Perfect adaptation means that the steady state response of the element is 
independent of the signal strength. Such a behaviour is typical of chemotactic 
systems (including human sense of smell) which initally responds to an abrupt 
change in interacting attractants and then settles down into a steady response 
[Painter K.J and Hillen T.(2002)], [Levchenko A, Iglesias P(2002)] have used 
a mechanism of this sort to model phosphoinosityl signaling in slime mold 
cells and neutrophils. The equations representing such a system are 



"' A: k l S-k 2 X l X 2 (6) 



dt 
dX 2 



= k 3 S - k 4 X x 



dt 

k 1 = k 2 = 2,k 3 = k 4 = 1,S = 1.2 



In some situations, like Frog oocyte maturation in response to progesterone 
and Apoptosis the switching is one way. At a critical signal strength we get 
an irreversible transition to large response starting from small responses. The 
transition can be activating or inhibitory giving rise to mutual activation or 
inhibition 



Mutual Activation Switch 



The equations representing such a system are 



dX 

—L = k E(X 1 )-k 2 E(X 1 )X 1 (7) 
dt 

E(X 1 ) = GBK(k 3 X 1 ,k 4 ,J 3 ,J 4 ) 

k = 0.4, fci = 0.1, k 2 = k 3 = 1, k 4 = 0.2, 
J 3 = J 4 = 0.05, S = 9.0 



The Goldbeter- Koshland function [Goldbeter A, Koshland DE(1981)] is graded 
and reversible. By graded we mean that the response increases continuously 
with signal strength. A slightly stronger signal gives a slightly stronger re- 
sponse. Reversible means that if we go from an initial value to a final value. 
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Although continuous and reversible, a sigmoidal response is abrupt. Tyson 
[Tyson (2003)] compares this to buzzer or a laser pointer, to activate the re- 
sponse one must push hard enough on the button, and to sustain the response 
one must keep pushing. When one lets up on the button, the response switches 
off at precisely the same signal strength at which it switched on. 

However, even in the presence of such a reversible term, the above switch 
shows irreversible behaviour due to a feedback loop that is present. 

Mutual Inhibition Switch 

In this switch, if signal is decreased enough, the switch will go back to the 
off-state. For intermediate stimulus strengths the response of the system can 
be either small or large, depending on how S was changed. This switch shows 
hystersis and is found in the lac operon in bacteria 

[Laurent M and Kellershohn N (1999)], the activation of M-phase-promoting 
factor (MPF) in frog egg extracts [Novak B and Tyson J(1993)], the autocat- 
alytic conversion of normal prion protein to its pathogenic form 
[Kellershohn N and Laurent M (2001)] and budding yeast cell cycle. This be- 
haviour has been observed in a number of experiments The equations repre- 
senting this system are 



-^ = k + hS- k 2 X l - k' 2 E(X 1 )X 1 (8) 

E(X 1 ) = GBK(k 3 ,hiX 1 ,J 3 ,Ji) 

k = 0.0, ki = 0.5, k 2 = 0.1, k 3 = l,k 4 = 0.2, 
J 3 = J 4 = 0.05, k' 2 = 0.5, S = 10.0 



Negative FeedBack Switch This type of regulation, commonly employed in 
biosynthetic pathways, is called homeostasis. 

The equations representing this system are 



dX 

- 1 ± = k E(X 1 )-k 2 SX 1 (9) 

E(X 1 ) = GBK(k 3: hX 1 ,J 3 ,J 4 ) 

k = 1.0, k 2 = 1.0, k 3 = 0.5, k 4 = 1.0, 
J 3 = J 4 = 0.01,5 = 0.1 



Activator Inhibitor Oscillator 

The classic example of an activator-inhibitor system is cyclic AMP production 
in the slime mold, Dictyostelium discoideum [Martiel J and Goldbeter A(1987)]. 
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External cAMP binds to a surface receptor, which stimulates adenylate cy- 
clase to produce and excrete more cAMP. At the same time, cAMP-binding 
pushes the receptor into an inactive form. After cAMP falls off, the inactive 
form slowly recovers its ability to bind cAMP and stimulate adenylate cyclase 
again. This mechanism lies behind all the curious properties of the cAMP 
signaling system in Dictyostelium: oscillations, relay, adaptation, and wave 
propagation. 

The equations representing this system are 



dX 

— - = k E(X 1 ) + kiS - k 2 X x - k' 2 X 1 X 2 (10) 
at 

— r— = k 5 Xi — k 6 X 2 
cit 

E(X 1 ) = GBK(k 3 X 1 ,k 4 ,J 3 ,J 4 ) 

k = 4.0,h = k 2 = k' 2 = k 3 = k 4 = 1.0, h = 0.1, k 6 = 0.075, 
J 3 = J 4 = 0.3, S = 0.2 



Substrate Depletion Oscillator 

This describes MPF oscillations in frog egg extracts [Novak B and Tyson J(1993)]. 
MPF is a dimer of a kinase sub-unit, cyclin-dependent kinase 1 (Cdkl), and a 
regulatory subunit, cyclin B. As cyclin B accumulates in the extract, it com- 
bines rapidly with Cdkl (in excess). The dimer is immediately inactivated by 
phosphorylation of the kinase subunit X can be converted into active MPF by 
a phosphatase called Cdc25. Active MPF activates Cdc25 by phosphorylating 
it. 

The equations representing this system are 
dX 

—± = k 1 S-(k' + k E(X 2 )X 1 (11) 
at 

dX 

-^ = (k' + hEiX^X, - k 2 X x 

E(X 1 ) = GBK(k 3 X 1 ,k 4 ,J ;i ,J 4 ) 

k' = 0.01, k = 0.4, k l = k 2 = k 3 = 1.0 
k 4 = 0.3, J 3 = J 4 = 0.05, S = 0.3 

where the GBK(u,v, J,K) is the Goldbeter-Koshland term 

GBK(u, v, J, K) = 2uk (12) 

v — u + vJ + uK + Jv — u + vJ + uK — 4{v — u)uK 
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Thus we can model a number of situations by coupling together these motifs. 
In the current work, we choose a pair of any of these elements to represent 
the types of nodes in our network. The parameters specified are fixed at these 
values following [Tyson (2003)]. The nodes are evolved by integrating the equa- 
tions using a fourth order Runge-Kutta routine for 100 iterates with time step 
0.1 for the switches and 0.001 for oscillators. 

2.3 Evolving the Network 

We follow the algorithm outlined in flow chart of Fig. 1. The algorithm evolves 
each node for a fixed time, updates the number of its neighbours Eq. (3), Eq. 
(4), update the current state by coupling it to the neighbours by Eq. (1), 
obtain the current position of the node by Eq. (2) and finally update the 
network by splitting the in and out degree. This consitutes one iteration of 
the network dynamics. We set the simulation as follows: we choose two of the 
above six types of nodes and set them up according to an initial golden mean 
network. The coupling coefficient e for this particular simulation is considered 
same for every node, though in principle the equations allow us to have a 
different value of e for every node. The offspring of each node is of the same 
as parent. We choose all possible pairs of the node to do the simulation. 



3 Characterisation and Analysis 

To understand the behaviour of a network a number of characterizing quanti- 
ties have been introduced in recent times 

[Erdos et al(1961); Barabasi(2001); Watts(1998); Newman]. It is important to 
choose the set of characterisers which reflect the essential features of the net- 
work with respect to the problems of interest. We classify the characterisers 
as the following : Topological characterisers which contain the information 
about the underlying topological structure of the connections in the net- 
work, Geometric characterisers , the growth or Entropic characterisers and 
characterisers for the dynamics. These characterisers are listed in Table 1 for 
three representative sets of parameters. 

3.1 Topological characterisers 

Degree Distributions 

The degree is an important characteristic of a node. Based on the degree 
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of the nodes, it is possible to construct measurements for the network. One 
of the simplest is the maximum degree which is simply the maximum con- 
nections that any node in the network has. The degree distribution, P(k) 
which expresses the fraction of nodes in a network with degree k contains 
more information. One may be interested in finding out if there is a corre- 
lation between the degrees of different nodes. Such correlations were found 
to have an important role in many network structural and dynamical prop- 
erties [Erdos et al(1961); Barabasi(2001); Watts(1998); Newman]. The most 
common choice is to find correlations between two nodes connected by a link. 
This correlation can be expressed by the joint degree distribution P(k, k'), i.e., 
as the probability of a link connecting two nodes of degree k and k! . Another 
way to express this is by giving the conditional probability that an arbitrary 
neighbor of a node of degree k has degree k' 

P{k\\k') = (13) 



This distribution gives a very detailed description of node degree correlations 
but, for fat tailed degree distributions as in scale- free networks, it is diffcult to 
evaluate experimentally, because of the poor statistics. A measure with better 
statistics is to computing the mean degree of the neighbors of nodes with a 
given degree, given by 

k nn {k) = Y, k ' P (k'\k) (14) 

k' 



In the simulation we use Eq. (14) for calulating vertex degree distribution. The 
vertex degree distribution can be calculated for the self degrees and neighbour 
degrees separately. For a random or small worlds network, the vertex degree 
distribution is a Poisson distribution in a large N limit with a peak at a 
particular degree bein most prominent. In contrast, for the scale free network 
of Barabasi and Albert, the degree distribution follows a power law for large 
k. In the network model that we consider we can get both types of degree 
distribution depending on parameter regimes. This is demonstrated in Figs 
2(a) and 2(b). the Fig 2(a) shows a typical small worlds behaviour whereas 
2(b) shows a power law like decay in large k region. Thus we can model 
both kinds of phenomena with the same model. As one changes the coupling, 
we get a transition from small worlds phase to a scale free phase. Further 
investigations of this transition are in progress. 

Clustering Coefficient 



A characteristic of the Erdos-Renyi model is that the local structure of the 
network near a node is a tree. More precisely, the probability of loops involving 
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a small number of nodes tends to in the large network size limit. This is in 
marked contrast with the profusion of short loops which shows up in many real- 
world networks. One way to characterise the presence of such loops is through 
the clustering coefficient. Two different clustering coefficients are frequently 
used, we can define a coefficient 

C ^ 

Nt — dijdjkdik 

i,j,k 

N 3 =Y1 a ij a jk 

(15) 



This coefficient gives the ratio of number of triangles i.e. the triples which are 
linked to each other to number of connected triples (i.e. three nodes connected 
to each other directly or indirectly). If we denote by Zj the number of links 
between neighbours of node i and hi is the number of neighbours of node i 
then we can also write clustering coefficient as 



a 
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ki(ki 1) 
N 4" 



(16) 



The difference between the two definitions is that Eq. (15) gives the same 
weight to each triangle in the network, while Eq. (16) gives the same weight 
to each node, resulting in different values because nodes of higher degree are 
possibly involved in a larger number of triangles than nodes of smaller degree. 
In the simulation result we calculate the coefficient given by Eq. (16). 

Centrality 



Greater the number of paths in which a node or link takes part, the greater the 
importance of this node or link for the network. Assuming that the interactions 
follow the shortest paths between two nodes, it is possible to quantify the 
importance of a node in this sense by its betweenness centrality 

(17) 



where a(j, i, k) is the number of shortest paths between nodes j and k which 
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pass through node i and o~(j, k) is total number of paths between nodes j and 
k. The sum is over all distinct pairs (j,k). Other centrality measures are given 
by closeness centrality D , network centrality N and stress centrality S 



A 




N 



S 



(max j d^) 



(18) 



where d^ = shortest distance between i,j 
3.2 Geometric characterisers 
Geodesic Lengths 

In the general case, two nodes of a complex network are not adjacent. In fact, 
most of the networks of interest are sparse, in the sense that only a small 
fraction of all possible links are present. Nevertheless, two non-adjacent nodes 
% and j can be connected through a sequence of m links (i, kl), (kl, k2), . . 
. , (km-1, j); such set of links is called a path between i and j, and m is the 
length of the path. We say that two nodes are connected if there is at least 
one path connecting them. 

A geodesic path (or shortest path) between nodes i and j, is one of the paths 
connecting these nodes with minimum length. The length of the geodesic paths 
is the geodesic distance d^ between nodes i and j. There can be more than 
one geodesic paths between two nodes. We take the mean geodesic distance 
as a characteriser of the network. The definition has a problem if the network 
contains unconnected nodes therefore this measurement is only performed 
after cleanup operation. The Dijkstra algorithm is used for calculating the 
shortest path on the network. 

Network Efficiency 

Network efficiency measure the ease of communication on a network. This can 
be defined as 



This quantity is also useful in characterizing robustness of network to attack 



i 



£4 



-i 



(19) 



N(N-l) 
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(see below) 



3.3 Entropic characterisers 
Growth rate and Entropy 

The entropy of a network can be calculated by using the following result 
based on Frobenius - Perron theory which we state without proof, for proof 
see [Lind D. and Marcus B. (1995)] 

Theorem : 

// G is an irreducible graph then h(Xc) = log\A{G) where h(Xa) is the entropy 
of the shift dynamical system X G associated with the graph Xa(g) ^ the Perron 
eigenvalue of the graph( i.e. in this case the largest real eigenvalue of the 
adjacency matrix). 

In general, however the above result provides an upper bound to entropy and 
we calculate this by diagonalising the adjacency matrix and finding the largest 
eigen value. Since our network is an evolving hypergraph , we compute this 
quantity at every time step and plot it with time in Fig 3. We find that 
the entropic bound saturates to certain fixed value (see Table 1) for a set 
of parameters and grows without bound for another set of parameters. This 
means that the network attains certain optimal size for a set of parameter 
values and has unrestrictive growth for another set. Continuously varying the 
parameters we have a transition from restrictive to unrestricted growth. 

Generalised Entropy 

Our network starts with nodes having a degree <=2. the degrees are added by 
splitting and by merger and removed by nodes dying. In the parameter regime 
where we get scale free behaviour the degree distributions can be fitted by 

P(> k) = et 2)/K 

e x q = [l + (l-q c )x]^ 
for 

1 + (1 - q c )x > (20) 

where K is related to the the value q c which best fits the data. This quantity is 
determined as follows: we take the logarithm of the probability for a series of 
different values of q c ranging from -3 to 3 in steps of 0.1. The values obtained 
from the q -exponential form are then fitted to data obtained from the model 
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by a linear fit. The value of q which gives the best fit is determined as value 

Qc 



This form is taken from the formulation of a non-extensive genrealised entropic 
function introduced by Tsallis [Gell-Mann M. and Tsallis C.(2004) ] 

The entropic function 



extremising this as described in [Gell-Mann M. and Tsallis C.(2004) ] gives 



this equation is identical to Eq.(20) if we assume q c = 1/(2 — q) and K = 
1/(2 — q) (3. Thus this non-extensive entropic function reproduces the behaviour 
of probability distribution. Therefore, it is probable that the correct entropic 
form for these models is non-extensive Tsallis entropy rather than by Shan- 
non - Boltzmann -Gibbs entropy. The consequences of this conjecture require 
further investigations. 



4 Response to Attack 

One of the key features of any network is its ability to withstand various kinds 
of attacks. An attack or damage to a network typically hampers its perfor- 
mance which leads to a decrease in its efficiency of communication or transfer 
of signals along various available paths in the network. Self organised net- 
works have various responses depending on the types of attack. Needless to 
say, this characteristic is of utmost importance for biological contexts where 
the response to an external attack could determine the survival, recovery or 
viability of a system. A number of attack strategies are usually considered 
while evaluating the response. These include targeted attacks (where a par- 
ticular node or connected segment is attacked based on some of its properties 
like vertex degree or length), multiple attacks (where random nodes are tar- 
geted), isolated attack (where a single node is targeted) etc. In the current 
simulation we consider isolated and multiple attacks. The measure of response 
to attack is given by network efficiency as defined in Eq. (19). The response of 
the network to these attacks is shown in the table where the average change 
in network efficiency with the number of nodes attacked is shown. This quan- 
tity is computed by choosing many nodes either one at a time or in a bunch 




l-^dk[p{k)f 
l-q 



(21) 




(2-9) 



(22) 
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and averaging over the decrease in network efficiency. As can be seen from 
the table, the network is more susceptible to multiple random attacks rather 
than to isolated single ones. This is a further manifestation of the fact that 
the spatio-temporal patterns found in the network are robust with respect to 
external perturbation such that removal of an element from a local pattern 
does not impair the network as much as removal of multiple randomly placed 
elements does. This is also indicative of the fact that there are more than one 
signalling paths available in any local neighbourhood and unless one blocks a 
large number of them the efficiency is not significantly impaired. 



5 Conclusion 

Therefore, in this paper we have proposed a generic model reproducing qual- 
itative behaviour of collective biological systems. The network incorporates a 
new feature of mobility which has been hitherto missing from such models. 
This property is shown to be essential in reproducing behaviour in presence of 
diffusion which many contexts show. The network goes beyond the assump- 
tion of boolean elements. It is shown that some of the shortcomings found in 
boolean networks can be overcome by making the state of an element a con- 
tinous variable. Prominient among such enhancements is an ability to model 
highly heterogenous and diverse environment where the internal state deter- 
mines and is determined by its surroundings, as in a cellular interaction mech- 
anism. At the same time the network allows us to have a number of states 
with switching attractors. Thus this model mimics the "robust yet fragile" 
nature of biological systems. 

This explicit coupling between intra and inter dynamics also gives an in- 
sight into several topological features of such networks. In particular, we have 
demonstrated that starting from a particularly simple initial configuration and 
without explicitly assuming any arbitrarily constructed rule for adding nodes, 
we can obtain both scale free nature as well as networks containing impor- 
tant hubs. The latter being reminiscent of small world networks. We present 
evidence for a transition between two regimes upon variation of parameters. 

We also show that interesting cooperative phenomena like community struc- 
tures and synchronisation emerge in these types of networks. This is clearly 
seen in Fig. 4 showing the spatial profile of internal variables. In this figure the 
regions of closely linked variations are distinctly separate from other regions. 

To investigate robustness of such a network to external perturbations we 
present a study of two different types of attack on the system. It is found 
that the multiple random attack is more effective than an isolated one. In 
current networks the reponse to the attack is reflected by a decrease in its 
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efficiency. 



The growth and size of the network follows an entropic law which shows sig- 
natures of non-extensivity. This leads to the speculation that Tsallis statistics 
is probably more relevant to biological networks. 

Several of these interesting features need further exploration and investigation 
which are currently under progress. However, in the present paper we only 
report that many qualitative aspects of the diverse behaviour of cooperation in 
biological systems can be modeled by a single generic network model. Various 
aspects of this model can give us an insight into important questions like 
quorum sensing in bacteria, synchronisation of biological elements, emergence 
of heterogeniety and morphogeneis, tissue and cell growth, tumour growth, 
protein interaction pathways and biochemial networks. 

For the network theory, this paper presents a network which shows number 
of different, well known network characteristics in a single model. Therefore, 
this effort may help in unifying some of different pieces of knowledge which 
we hope would be useful in understanding general networks better. 

Both authors contributed equally to the work. Supplementary materials like 
programs, data set and additional figures are freely available from the authors 
upon request. 
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6 Figure Captions 

(1) Figure 1: Figure 1 describes the flow of the network evolution. The 
blocks in the figure describe how various stages of dynamics are handled 
and how they depend on each other. 

(2) Figure 2: Figure 2 describes the typical vertex degree distributions ob- 
tained. The figures 2(a) and 2(b) are generated by initially setting up a 
golden mean network with two different kinds of node,namely the Mutual 
activation switch Eq.(8) and the activator - inhibitor oscillator Eq.(ll). 
The network was evolved with value of R = 0.35, and the same values of 
node parameters as in [Tyson (2003)]. The value of the coupling strength 
e was 0.66 for 2(a) and 0.033 for 2(b). Each node was evolve with respect 
to internal dynamics by using a fourth order unge-Kutta method. The 
step size for the switch was 0.1 and for the oscillator 0.001. 

(3) Figure 3: Describes the entropic bound for growth as described in the 
paper. Fig. 3(a) describes the behaviour where the growth of the network 
saturates after some time and Fig. 3(b) where it grows without bounds. 

(4) Figure 4: shows the patterns found in internal profile of the network at 
various times (a)-(e). The presence of correlated patterns can be easily 
seen from this figure. 
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Set up a Golden Mean network and initialize 
the network and node parameters 

< 

I 

Evolve each node by integration using fourth 
order Runge-Kutta routine 

I 

Update the coupling between each node and 
its neighbors 

i 

Update the position and neighbors of each 
node 

i 

Evolve the network by state splitting 




* 

Calculate the network efficiency, vertex degree 
distribution, entropy and fractal dimension 



Fig. 1. 
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